# ###### ALTERNATIVE MECHANISM TESTS: ST PROPORTION ######

interact_theme = theme_set(
  theme_classic() +
    theme(axis.text = element_text(face = "bold", size = 16),
          axis.title = element_text(face = "bold", size = 18),
          legend.title = element_text(face = "bold", size = 18),
          legend.text = element_text(face = "bold", size = 16))
)

prop.st.covlabs = c(
  "ST Reservation", "Female Reservation", "Prop. ST",
  "Fem $\\times$ ST Reservation","ST Res $\\times$ Prop. ST",
  "Female Res $\\times$ Prop. ST", "Fem $\\times$ ST Res. $\\times$ Prop. ST")

controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

form = paste("caste_interact_bi_comp ~ vstreserv*vwreserv*prop_ST91 + " , controls)

# Run interaction model
m <- form %>% as.formula %>% lm_robust(data = reds_p3_rand, fixed_effects = ~districtid,
                                       clusters = villageid, se_type = "stata")
vcov1 = m$vcov
# Interaction values with and without women's reservations
no_vst_propST = reds_p3_rand$prop_ST91[reds_p3_rand$vstreserv==0]
vst_propST = reds_p3_rand$prop_ST91[reds_p3_rand$vstreserv==1]

# Calculate marginal effects
ame_no_vst = m$coefficients["vwreserv"] + m$coefficients["vwreserv:prop_ST91"]*no_vst_propST
ame_vst = m$coefficients["vwreserv"] + m$coefficients["vstreserv:vwreserv"] +
  m$coefficients["vwreserv:prop_ST91"]*vst_propST + m$coefficients["vstreserv:vwreserv:prop_ST91"]*vst_propST

# Calculate var-covar
no_vst_vcov = vcov1["vwreserv", "vwreserv"] + vcov1["vwreserv:prop_ST91", "vwreserv:prop_ST91"]*((no_vst_propST)^2) +
  2*no_vst_propST*vcov1["vwreserv", "vwreserv:prop_ST91"]


vst_vcov = vcov1["vwreserv", "vwreserv"] + vcov1["vstreserv:vwreserv", "vstreserv:vwreserv"] +
  vcov1["vwreserv:prop_ST91", "vwreserv:prop_ST91"]*((vst_propST)^2) +
  ((vst_propST)^2)*vcov1["vstreserv:vwreserv:prop_ST91", "vstreserv:vwreserv:prop_ST91"]+
  2*vcov1["vwreserv", "vstreserv:vwreserv"] +
  2*vst_propST*vcov1["vwreserv", "vwreserv:prop_ST91"] +
  2*vst_propST*vcov1["vwreserv", "vstreserv:vwreserv:prop_ST91"] +
  2*vst_propST*vcov1["vstreserv:vwreserv", "vwreserv:prop_ST91"] +
  2*vst_propST*vcov1["vstreserv:vwreserv", "vstreserv:vwreserv:prop_ST91"] +
  2*vst_propST*vst_propST*vcov1["vwreserv:prop_ST91", "vstreserv:vwreserv:prop_ST91"]
no_vst_se = no_vst_vcov %>% sqrt
vst_se = vst_vcov %>% sqrt

plot.dta = cbind.data.frame(prop_ST = no_vst_propST, ame = ame_no_vst,
                            vstreserv = 0, se = no_vst_se) %>%
  bind_rows(cbind.data.frame(prop_ST = vst_propST, ame = ame_vst,
                             vstreserv = 1, se = vst_se)) %>%
  mutate(ci.lo = ame-1.96*se, ci.hi = ame+1.96*se)

# Plot marginal effects
pl1 = ggplot(data = plot.dta, aes(x = prop_ST, y = ame), group = as.factro(vstreserv))+
  geom_line(aes(color = as.factor(vstreserv)))+
  geom_ribbon(aes(ymin = ci.lo, ymax = ci.hi, fill = as.factor(vstreserv)), alpha = 0.3)+
  geom_rug(aes(color = as.factor(vstreserv)), sides = "b", alpha = 0.3)+
  geom_hline(yintercept = 0, color = "black", linetype = "longdash")+
  scale_color_manual(name = "Caste Quotas", values = c("gray", "firebrick4"),
                     label = c("No", "Yes"))+
  scale_fill_manual(name = "Caste Quotas", values = c("gray", "firebrick4"),
                    label = c("No", "Yes"))+
  labs(x = "Proportion Scheduled Tribes", y = "Average Marginal Effects")+
  theme_pubclean()+
  theme(axis.title = element_text(size = 16),
        axis.text = element_text(size = 16),
        legend.text = element_text(size = 16),
        legend.title = element_text(size = 16))
pl1
ggsave(pl1, filename = paste0(fig.out, "FigureL12b.png"),
       width = 10.4, height = 7.6)

  
  
screenreg(l = list(m), include.ci = F, include.rmse = F,
            include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
            digits = 3, include.nclust = F,
            omit.coef = c("Intercept|^prop_ST91$|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12"),
            reorder.coef = c(2, 1, 3, 5, 4, 6),
            custom.coef.names = c("Women's Quota", "Caste Quota",
                                  "W X C Quota", 
                                  "W Quota X ST Mobilization",
                                  "C Quota X ST Mobilization",
                                  "W X C Quota X ST Mobilization"),
            stars = c(0.01, 0.05, 0.1),
            custom.model.names = c("REDS"),
            custom.gof.rows = list(
              "Mechanism Type" = c("ST Proportion"),
              "Caste Quota Type" = c("ST Quota"),
              "District FE" = c("Yes"),
              "Controls" = c("Yes")), 
            file = paste0(tab.out, "TableL25_Col1.tex"))  

# ###### ALTERNATIVE MECHANISM TESTS: ST PARTY VOTE SHARE ######
stres.scmob.covlabs = c(
  "ST Reservation", "Female Reservation", "ST Party Vote Share",
  "Fem $\\times$ ST Reservation","ST Res $\\times$ ST Party VS",
  "Female Res $\\times$ ST Party VS",
  "Fem $\\times$ ST Res. $\\times$ ST Party VS")

form = paste("caste_interact_bi_comp ~ vstreserv*vwreserv*SC_party_pc_vids + " , controls)

# Run interaction model
m <- form %>% as.formula %>% lm_robust(data = reds_p3_rand, fixed_effects = ~districtid,
                                       clusters = villageid, se_type = "stata")
vcov1 = m$vcov

# Interaction values with and without women's reservations
no_vst_propST = reds_p3_rand$SC_party_pc_vids[reds_p3_rand$vstreserv==0]
vst_propST = reds_p3_rand$SC_party_pc_vids[reds_p3_rand$vstreserv==1]

# Calculate marginal effects
ame_no_vst = m$coefficients["vwreserv"] + m$coefficients["vwreserv:SC_party_pc_vids"]*no_vst_propST
ame_vst = m$coefficients["vwreserv"] + m$coefficients["vstreserv:vwreserv"] +
  m$coefficients["vwreserv:SC_party_pc_vids"]*vst_propST + m$coefficients["vstreserv:vwreserv:SC_party_pc_vids"]*vst_propST

# Calculate var-covar
no_vst_vcov = vcov1["vwreserv", "vwreserv"] + vcov1["vwreserv:SC_party_pc_vids", "vwreserv:SC_party_pc_vids"]*((no_vst_propST)^2) +
  2*no_vst_propST*vcov1["vwreserv", "vwreserv:SC_party_pc_vids"]


vst_vcov = vcov1["vwreserv", "vwreserv"] + vcov1["vstreserv:vwreserv", "vstreserv:vwreserv"] +
  vcov1["vwreserv:SC_party_pc_vids", "vwreserv:SC_party_pc_vids"]*((vst_propST)^2) +
  ((vst_propST)^2)*vcov1["vstreserv:vwreserv:SC_party_pc_vids", "vstreserv:vwreserv:SC_party_pc_vids"]+
  2*vcov1["vwreserv", "vstreserv:vwreserv"] +
  2*vst_propST*vcov1["vwreserv", "vwreserv:SC_party_pc_vids"] +
  2*vst_propST*vcov1["vwreserv", "vstreserv:vwreserv:SC_party_pc_vids"] +
  2*vst_propST*vcov1["vstreserv:vwreserv", "vwreserv:SC_party_pc_vids"] +
  2*vst_propST*vcov1["vstreserv:vwreserv", "vstreserv:vwreserv:SC_party_pc_vids"] +
  2*vst_propST*vst_propST*vcov1["vwreserv:SC_party_pc_vids", "vstreserv:vwreserv:SC_party_pc_vids"]
no_vst_se = no_vst_vcov %>% sqrt
vst_se = vst_vcov %>% sqrt

plot.dta = cbind.data.frame(prop_ST = no_vst_propST, ame = ame_no_vst,
                            vstreserv = 0, se = no_vst_se) %>%
  bind_rows(cbind.data.frame(prop_ST = vst_propST, ame = ame_vst,
                             vstreserv = 1, se = vst_se)) %>%
  mutate(ci.lo = ame-1.96*se, ci.hi = ame+1.96*se)

# Plot marginal effects
pl1 = ggplot(data = plot.dta, aes(x = prop_ST, y = ame), group = as.factro(vstreserv))+
  geom_line(aes(color = as.factor(vstreserv)))+
  geom_ribbon(aes(ymin = ci.lo, ymax = ci.hi, fill = as.factor(vstreserv)), alpha = 0.3)+
  geom_rug(aes(color = as.factor(vstreserv)), sides = "b", alpha = 0.3)+
  geom_hline(yintercept = 0, color = "black", linetype = "longdash")+
  scale_color_manual(name = "Caste Quota", values = c("gray", "firebrick4"),
                     label = c("No", "Yes"))+
  scale_fill_manual(name = "Caste Quota", values = c("gray", "firebrick4"),
                    label = c("No", "Yes"))+
  labs(x = "ST-mobilizing Party Vote Share (Vidhan Sabha)", y = "Average Marginal Effects")+
  theme_pubclean()+
  theme(axis.title = element_text(size = 16),
        axis.text = element_text(size = 16),
        legend.text = element_text(size = 16),
        legend.title = element_text(size = 16))


pl1
ggsave(pl1, filename = paste0(fig.out, "FigureL12a.png"),
       width = 10.4, height = 7.6)

screenreg(l = list(m), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|^prop_ST91$|^SC_party_pc_vids$|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12"),
          reorder.coef = c(2, 1, 3, 5, 4, 6),
          custom.coef.names = c("Women's Quota", "Caste Quota",
                                "W X C Quota", 
                                "W Quota X ST Mobilization",
                                "C Quota X ST Mobilization",
                                "W X C Quota X ST Mobilization"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("REDS"),
          custom.gof.rows = list(
            "Mechanism Type" = c("ST Party Vote Share"),
            "Caste Quota Type" = c("ST Quota"),
            "District FE" = c("Yes"),
            "Controls" = c("Yes")), 
          file = paste0(tab.out, "TableL25_Col2.tex")) 



